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ABSTRACT 

Investigation  of  two-dimensional  transonic  flows  is 
extended  to  axi- symmetric  problems.  This  is  of  considerable 
practical  interest,  for  example,  with  regard  to  missiles  or 
aircraft  engines  which  approximate  much  more  closely  bodies  of 
revolution  than  two-dimensional  bodies.  The  main  concern  with 
axi -symmetric  flows  lies  not  only  with  the  complexity  of  the 
governing  nonlinear  partial  differential  equation  which  is 
mixed  of  elliptic-hyperbolic  type  but  also  with  the  lack  of  a 
general  method  for  accurately  solving  this  type  of  equation. 
We  solve  the  nonlinear  transonic  equation  using  separation 
of  variables  technique,  which  yields  two  nonlinear  ordinary 
differential  equations.  The  x- dependence  can  be  integrated 
numerically,  and  the  solution  for  the  r-dependence  can  be 
obtained  using  the  expansion  method  originated  by  Van  Dyke. 
This  works  well  with  only  three  terms  in  the  expansion.  The 
sonic  solution  of  these  equations  is  obtained  analytically 
since  both  equations  are  simple  enough  to  be  integrated  for 
this  case  (1^=1.0).  The  small  parameter  (l-M^2)  plays  an 
important  role  in  specifying  the  shape  of  the  boundary 
surfaces  for  external  axi -symmetric  steady  flow  of  interest 
for  design.  A  Navier- Stokes  solver  was  used  to  compute  the 
inviscid  flow  to  confirm  our  results  in  the  region  over  the 
surface  where  the  small  perturbations  apply. 
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I .     INTRODUCTION 

Transonic  flow  is  a  classical  problem  in  gas  dynamics  and 
yet  not  many  exact  solutions  are  available.  The  main 
difficulty  lies,  of  course,  in  the  nonlinearity  of  the 
equations  and  in  the  boundary  conditions . [Ref .  1] 

It  is  of  interest  to  extend  two-dimensional  flow  solutions 
to  axi- symmetric  problems,  which  are  more  practical  for 
designing  bodies  of  revolution  such  as  missiles  or  aircraft 
engines.  The  goal  is  to  design  a  shock- free  body  or  surface  so 
we  can  minimize  the  penalty  of  the  wave  drag  that  is  caused  by 
shock  waves.  The  small  perturbation  technique  may  be 
implemented  for  simplifying  the  governing  equations  and  the 
mathematical  procedures  for  solving  the  resulting  equations  in 
nonlinear  flows.  This  research  shows  a  solution  to  nonlinear 
two-dimensional,  external  flow,  over  an  afterbody  axi- 
symmetric  surface.  Beginning  with  the  transonic  equation  for 
an  axi -symmetric  body  in  Chapter  II  and  applying  separation  of 
variables  gives  two  nonlinear  ordinary  differential  equations 
which  lead  into  two  exact  solutions.  One  of  the  nonlinear 
ordinary  differential  equation  can  be  integrated  numerically, 
and  other  one  can  be  solved  using  the  outer  expansion  method 
by  Van  Dyke  [Ref.  2].  Chapter  III  shows  the  derivation  of 
pressure  coefficient  and  local  Mach  number  for  different 


surfaces.     Supersonic   boundary   surfaces   are   obtained 
implicitly  in  Chapter  IV  for  W0B=1 .  05  , 1 . 1,  and  1.2  and  they  are 

examined  by  a  3D  Navier- Stokes  solver,  run  as  an  Euler  solver 
for  this  axi- symmetric  problem  in  Computational  Fluid  Dynamics 
(CFD) .  In  addition,  the  sonic  surface  is  derived  analytically 
in  Chapter  V.  This  sonic  boundary  surface  is  also  examined 
using  3D  Navier- Stokes  solver  in  CFD  to  verify  the  results. 


II.     TRANSONIC  EQUATION  FOR  AXI- SYMMETRIC  BODY 

The   small   perturbation,   non-linear,   axi- symmetric, 
transonic  equation  is  written  as  [Ref.  3] 

U-l£)9„+±-gjUiz)   =<PX<PXX  (1) 


Where  the  non-dimensional  axial,  x,  and  radial,  r,  velocity 
potentials  have  been  defined  as  [Ref.  1] 

<px  =  bti    (y+1)-^,    (pr  =  Ml    (Y+D-^         <2) 


where  M^  is  free  stream  Mach  number  and  y  is  ratio  of  heat 
capacities. 

An  exact  solution  to  Equation  (1) ,   the  transonic  equation, 
can  be  found  from 

<p(x,r)  =  cps(x,r)  +  (l-Mi)x  (3) 

Which  patterned  after  the  two  dimensional  case  given  by 
Biblarz  [Ref.  1]  .  Here  <ps  (x,r)  is  a  separable  solution  for 
the  velocity  potential. 


The  transonic  equation  in  axi- symmetric  form,  Equation  (1) , 
may  be  separated  by  letting  the  potential  function  (p  (x,r) 
equal  to 


<p  (x,r)  =1  (x)  C(r)  +  (l--Mf)x 


(4) 


Substituting  the  above  function  into  the  transonic  equation 
results  in  two  ordinary,  second  order,  non- linear  differential 
equations 


dx   dx2 


(5) 


dr2       r  dr 


(6) 


where  X      is  the  separation  constant. 

The  solution  to  the  first  O.D.E. (5)  is  obtained  by  multiplying 

both  sides  by  d\/dx, 


dx 


dx  dx2  )      dx 


(7) 


or 


_d_ 
dx 


=  0 


(8) 


Thus 


dx  2     s 


Rearranging 


dx  = 


c£ 


(  i^2  +  a ) 


and  integrating 


(j  *«*♦«) 


(10) 


— r  (ID 


where  a  and  x0   are  integration  constants. 

The  solution  to  the  second  non-linear  O.D.E.(6),  is  obtained 
by  an  outer  expansion  method.  [Ref.  2] 

C(r)=-i-  +  (l-ndjfj  (r)  +  (l-MZ)2f2(r)   +..  .     (12) 
Ar2 

Where  (1-wf)  represents  a  small  parameter  and  the  first  term 
is  the  purely  sonic  solution. 


By  taking  first  and  second  derivatives  and  substituting 
them  into  the  original  2nd  order,  non-linear,  O.D.E.(6),  the 
solution  becomes 


C(r)  = 


kr: 


4  +  |l-M„2|aXr(^+2)+il^l!a2X2r2(^+2)  + 

28+8/8 


(13) 


where  a   is  constant 


Equations  11  and  13  are  solutions  to  the  differential 
equations  5  and  6.  Boundary  conditions  are  needed  to 
determine  the  constants  a  ,  C1 ,    a,    and  X  . 

The  constants  Cx ,    a,    and  X  can  be  shown  to  be  related  by  the 

expression  [Ref.l] 


(  C   N 
{   A 


1.2426 


0.7574 


=  1.08x10 


(14) 


and 


a  =  +  Cx  1 1  -Mt 


2  11.7574 


(15) 


The  positive  sign  above  is  for  AC  2  1.0  and  the  negative  sign 
is  for  M^,  <  1.0. 


Upon  inserting  Eqn. (15)  into  Eqn. (11) ,  we  have 


XX°~l 


dl 


1  XZ2  +  a  \±-mZ  I1-7574 


(16) 


Factoring  out      Cx ,    we   obtain 


X  x0  -  Cx 


1 


dt 


Ixli  +  Il-M*  11-7574 


2        C 


(17) 


Introducing  new  variables 


E  -  I 


3.X 

V2C, 


(18) 


^  (v/8+2)  =  r  (v/^  +  2)  ^ 


(19) 


Equations    (17)    and    (13)    become 


x    xQ  -  Cx 


W2  ( — - — i 

J  [52+|i-i^|1-7574]3 


(20) 


c  =  _(aX)0-4142 


X? 


P2 


r2\   2 


4+  1-AL  f  (v8+2)+ r2 

1     '         50.63 


Finally,  defining  two  new  variables 


(21) 


x~=C^   x{f^ 


(22) 


I   - 


C  X 


(a  X) 


0.4142 


(23) 


Equations    2  0   and   21   become 


J  [S^li-j^l1-7574]3 


where  x0=0 


(24) 


f  =  ^ 


f2 


4  +  ll-M2  I  I   <v^  +  2)  +_0_JO_  ~2(v^+2: 

1    "'  50.63 


(25) 


Equation  (24)  will  be  numerically  integrated  and  plotted  in 
Fig.  (1)  as  I  (x)  vs.  x  ,  and  Equation  (25)  will  be  evaluated 
and  plotted  in  Fig.  (2)  as  X,  (?)    vs.  ?    for   Af„=  1.05,  1.1, 

1.2.  These  results  will  be  used  later  to  obtain  transonic 
surfaces  with  their  corresponding  Cp    and  local  Mach  number 


profiles.  In  Fig.  (2)  ,  Equation  (25)  is  plotted  up  to  its 
minimum  value  with  increasing  f.  Thereafter,  a  constant 
value  is  patched  in  to  the  right  of  the  +  mark.  This  is  done 
because  of  the  need  to  keep  £  constant  reflecting  the 
necessary  behavior  of  the  function  for  range   f  [Ref .1] . 


III.    PRESSURE  COEFFICIENT  AND  LOCAL  MACH  NUMBER 

The   pressure   coefficient   for   a   flow   with   small 
perturbations  is  given  by  [Ref.  4] 


Cp=  - 


2u 


u. 


\2 


(26) 


The  linearized  pressure  approximation  for  axi- symmetric  flow 
is 


"2  ux 
r     =  - 


(27) 


Recall  the  axial  velocity  potential,  Equation  (2) ,  Chapter  II, 


u. 


<px  =  «C  (Y  +  D  -— 


(28) 


thus 


CL    ^2 


<P; 


at.:  (Y+D 


(29) 


10 


substitute  Equation  (29)  into  Equation  (27) ,  yields 


C  =    ~2  ** 
P   Ad  (Y+D 


(30) 


We  need  to  take  the  derivative  of  the  potential  function  with 
respect  to  x  in  equation  (4) . 


•*■ c  il +  (1  ~^2) 


(31) 


Rewriting  equation  (9)  with  the  constant  a  inserted  becomes 


dL    - 
dx 


3X 


52  +  q  |  (1-MZ) 


_M2\     11.7574 


,  1 


(32) 


Recall  Equation  (18) 


I   =  I 


2L) 


(33) 


Inserting  Eqn.(33)  into  Eqn. (32)  and  factoring  Cx   out,  yields 


rt    - 
dx 


=  cx2  [  r  +  |(i-m.2)  i1-7574  ] 


(34) 


11 


Recall  Equation  (23)  and  rearrange 


c  =  (a   X)0-4142  ^ 


(35) 


Rewriting  Equation  (31)  the  potential  function  as 


cpx=  (aX)x0'4142  £  cx3  [^Hd-^2)!1-7574]3^!-^2)   <36) 


Thus 


<P* 


a  0.4142/-.  3 


0.5858 


i-  f  [  E2  +  |(l-^)|1,7574]3  +  (l-^2) 


(37) 


Recall    Equation    (14)    Chapter   II 


/^  \ 


V  a  / 


1.2426 


0.7  57  4 


=  1.  08xl0-2 


(38) 


then   finally   Equation    (37)    becomes 


<px  =  0.2208  f  [  £2   +    Kl-M2)  I1'7574  ]3    +    (1-i^2) 


(39) 


Rewriting  the  coefficient  of  pressure 


C  = 


-2 


p  Ml(y+i) 


0.220S£[£2  +  |(1-M2)  I1-7574]3 +  (1-2^ 


(40) 


12 


The  previous  expression  for  Cp    will  be  used  further   in 

determining  the  local  Mach  number. 

The  perturbation  velocities  are  assumed  to  be  vary  small 
compared  to  the  undisturbed  uniform  velocity  £/_.  [Ref .  5] 

Hence  at  a  point  near  the  body,  the  velocity  vector  V   for 
2-D  flow  is  given  by 

v  =  i(um  +  ux)    +  j  uz  <41) 

We  can  obtain  the  local  Mach  number  M  in  terms  of  the 
perturbation  velocities  and  the  local  speed  of  sound  c  as 


\v\2    _     (Um   +  ux)2   +  ul  (42) 


AT  = 


C2  C2 


Neglecting  the  higher  order  ratios  of  the  free  stream  to 
perturbation  velocities,  since  they  are  considerably  smaller 
than  unity,  we  have 

„2  -  ui  (1+i£) 

c2 


13 


The    ratio    of    the    local    speed    of    sound    to    the    undisturbed 
uniform  speed  of   sound  becomes 


c2    _     T 
cl         Tm 


(44) 


or 


•2  T0/TK 


cl  To/T 


(45) 


Where    TQ    is   the   stagnation   temperature. 

We   can  write   temperature   ratios   as    functions    of   Mach  numbers 


as 


2       i  +  M^Lid 

—,  =  r^n —  (46) 

d       i  +    (Y-DM2 

2 


Substituting  Equation  (45)  into  Equation  (46)  and  rearrange, 
we  have 


'1  +  ^\1  +  ±LZ11M2 


M2   =  Ml 


2 


14 


from  Equation  (27) 

-2  u 


Cp  =   ^p  (48) 


Rearranging  and  solving  for  local  Mach  number 


M2     m  ^      (1-Cp) 


+  Jtzl 


wf  c 


(49) 


This  equation  will  be  used  to  show  the  local  Mach  number  for 
each  surface  at  Mm   =   1.05,  1.1  and  1.2  [Ref.  6]. 


15 


IV.     SUPERSONIC  BOUNDARY  SURFACES 

The  boundary  conditions  which  need  to  be  applied  require 
that  the  gradient  of  q>  vanish  far  ahead  of  the  body  and  that 
the  flow  be  tangential  to  the  surface  [  Ref .  7] .  In  terms  of 
perturbation  velocities  the  boundary  conditions  become 

(#)      =  -  (50) 

y^  I  surface  3. 

Recall  the  modified  velocity  perturbation  potential  Eqn. (2) 

cpr  =  Mi    (y+l)ik  (51) 


then,  we  have 


"  — = ~ (52) 

Ml    (Y+D 


by  substituting  Equation  (52)  into  Equation  (50) ,  we  obtain 

^   -  *'  <53, 


dxl       mI   (y  +  1) 


16 


By  differentiating  Equation  (4)  with  respect  to  r  gives 


<Pr  =  I 


dr 


(54) 


Recalling  Eqn. (18)  and  substituting  Eqn. (21)  into  Eqn. (23) , 
and  taking  the  derivative  with  respect  to  f,  the  above 
equation  becomes 


<pr  =  0.0848  I 


dr 


(55) 


Substituting  Equation  (55)  and  (13)  into  equation  (53),  we 
have 


dr 
dx) 


0.0848 


I 


-4    +  2.8284    ll-Afil   f1-8284 


ac(y+D      Lr 

+  0.1512    {1-mZ)2  f6>657 


(56) 


It  can  be  noticed  here  that  \   and  £   are  given  implicitly  as 

function  of  x   and  r  respectively. 

From  Eqns.(19)  and  (22),  we  can  rearrange  Equation  (56)  as 


dr         3  .26x10 
dx 


I 


-8 

S*3 


+  2.8284    \l-ld\   f 1-8284 


+  0.1512    (1-Mi)2  f6-657 


(57) 


17 


Further  arrangment  of  the  previous  equation  yields 


dr 


—  -2.8284|l-Ad|   f1-8284   -  0.1512  (1-Afc2)2  f6-657 


f3 


(58) 


-3  .26X1Q-2 


I  dx 


Now,  we  can  compute  the  exact  boundary  surfaces  out  of  the 
equation  above.  Starting  with  the  left  hand  side  of  the 
equation 


Z(f)=- 


f3 


8-2 .8284  1-AC  f 


.VfZ  I  ?4.8284_ 


0.1512 (l-AC)2f 


_M-2x  2^9.657 


(59) 


By  defining  equation  (59)  we  can  plot  Z  (f)  vs.  f  in  Fig. 
(3)  for  Mm   =  1.05,  1.1,  and  1.2. 

The  minimum  values  of  iQ    for  different  Mach  numbers  can  be 

expressed  as  fQ   =  r0  (Mj      which  has  been  found  to  be  [Ref .  1] 


ro   = 


1.2074 

1-^210.2071 


(60) 


18 


Fig.  (4)  represents  f0   vs.  M„   and  shows  a  symmetry  close  to 
MM=1.0  where  f0  =  °°.    So  for  our  purpose  to  represent  the 
condition  sufficiently  close  to   Af,.  =  1.0,  we  will  choose  a 
finite  value  of  f0    (perhaps  as  4.0). 
Defining  the  right  hand  side  of  Equation  (58)  as 


K(x)  = 


-3.26X1Q-2 
Ad(Y+D 


I  dx 


(61) 


taking  the  derivative  of  Equation  (24)  and  substitute  it  into 
Equation  (61) ,  yields 


K{x)  = 


-3 .26x1 


Ml{y 


k10~2  C    £d| 

+  1)    Jo   [J2+|l-^2|l-7"4]3 


(62) 


integrating  by  parts,  we  obtain 


K(x)    = 


-2  .45X10"2 


(y+dm: 


(^-li-Mii1-7574)3*!!-^2!1 


1716 


(63) 


the  equation  above  will  allow  us  to  plot  K(x)vs.x  in 
Fig. (5)  for  M„   =  1.05,  1.1,  and  1.2.   Where  K (x)  <  0 

for  MM2  1.0  and  K (x)  >   0  for  Mm   <1.0. 
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Therefore,  integrating  both  sides  of  equation  (58)  becomes 


I 


f3 


df = 

-.1512  (1-Mj)2r6-657 

(64) 


_§_-9      R9RA  ll  _\f2lf!l-8284_     i  m  9  fi  .u2\2f  6.657 


-2.82  84  1-AC  f1, 8284-.  1512  (1-M„)2f 


x 


-3.26X10-         ldjt 


Ad(Y+l) 


Equation (64)  is  numerically  integrated  and  which  then  allows 
us  to  determine  the  boundary  surfaces  in  dimensionalize  and 
non-dimensionalize  (normalize)  form  as  shown  in  Figs.  6  and  7 
for  Mm    =    1.05,  1.1,  and  1.2.  Now,  we  can  determine  Cp  from 

Eqn. (40)  and  local  Mach  number  from  Eqn. (49)  for  the  three 
surfaces  obtained  at  different  Mach  numbers. 

In  Figs.  8-13,  the  supersonic  boundary  surfaces  for 
1^=1.05,1.10,  and  1.20  are  plotted  with  their  pressure 
coefficients  and  local  Mach  profiles  in  dimensionalized  and 
non-dimensionalized  (normalized)  form.  These  results  will  be 
verified  later  using  a  3-D  Navier- Stokes  solver  in 
Computational  Fluid  Dynamics  (CFD) . 
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V.      SONIC  BOUNDARY  SURFACE 

The  sonic  flow  (Afoo=1.0)  can  be  derived  in  explicit  form 

from  Equations  (24)  and  (25) .  The  resulting  equations  for 
sonic  flow  are 


l"4  (65) 


/ 


and 


I   =  A:  (66) 


i2 


Rearranging  Equation  (65) 


J   o 


z    3  dz  (67) 


since  Equation  (67)  is  integrable,  we  obtain 


x  =   3  I 


or 


<?3 
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(68) 


I   =  ^-  (69) 
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We  can  write  Equations  (11)  and  (13)  in  sonic  form  as 


x 


I 


(1 » 1" 


(70) 


and 


C  -  ji-i  (71) 

A  r2 


integrating  Equation  (70) ,  we  have 


l(x)    =   -L  X  x3  (72) 

1 8 


From  Equation   (4)  ,   we  can  write  the  sonic  perturbation 
potential  as 


cp(x,r)  =  $U)C(r)  (73) 


Substituting  Equations  (71)  and  (72)  into  Equation  (73)  yields 

(p(x,r)  =  -|*1  (74) 

9  x2 


Equation  (57)  can  be  written  for  Mm=1.0    as 

di  -  0 .2608  I  (rjs) 

dx  (y  +  1)  f3 
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Rearranging  Equation  (75) ,  we  have 


f3 dr  =-0.1087  I  dx  <76) 


Therefore,  integrating  both  sides  of  Equation  (76)  ,  we  obtain 


|  r3  dr  =-0.1087  /  I  dx  (77) 


Inserting  Equation  (69)  into  Equation  (77)  yields 


f4  -f04  =  -0.0040X4  (78) 


where  f 0  =4 . 0  has  been  chosen  as  an  asymptotic  value  to 

represent  the  condition  sufficiently  well  at  Af^l . 0  .  Equation 

(78)  will  allow  us  to  determine  sonic  boundary  surface. 
The  pressure  coefficient  can  be  written  from  Eqn. (40)  as 


C    -   -0-4416>  r|  (79) 

p    (Y  +  D    * 


in  case  of  AL=1 . 0 . 
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The  local  Mach  number  can  be  written  from  Equation  (49)  for 
the  sonic  surface  as 


M     -, 733; — r  (80) 


The  sonic  surface  is  represented  dimensionally  in  Fig.  (14) 
with  the  pressure  coefficient  and  local  Mach  number  profiles, 
and  normalized  in  Fig.  (15)  with   f 0  =4 . 0  for  (M^l.O).  The 

sonic  boundary  surfaces,  in  both  forms,  are  plotted  with  the 
supersonic  boundary  surfaces  in  Figs.  16  and  17. 
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VI.     COMPUTATIONAL  FLUID  DYNAMICS  (CFD) 

The  history  of  Computational  Fluid  Dynamics  (CFD)  is 
closely  tied  to  the  rapid  advances  in  digital  computer  which 
has  a  great  impact  on  problems  of  design  in  modern  engineering 
practice.  These  problems  can  now  be  solved  at  very  little  cost 
in  a  few  seconds  of  computer  time  which  would  have  taken  years 
to  work  out  with  the  computational  methods  and  computers 
available  twenty  years  ago.  [Ref.  8] 

The  CFD  will  be  used  to  compute  the  axi- symmetric  flow 
over  the  boundary  surfaces  obtained  by  the  small  perturbation 
method. 

The  computer  program  GRAPE  [Ref. 9],  an  acronym  derived 
from  Grids  about  Airfoils  using  Poisson's  Equation,  has  been 
written  by  R.  Sorenson  at  Ames  Research  Center  to  generate 
two-dimensional  finite  difference  grids  about  airfoils  and 
other  shapes  by  the  use  of  the  Poisson  differential  equation. 
Outer  and  inner  boundaries  are  specified  for  the  C-type  grid, 
where  the  surface  of  the  body  is  treated  as  the  inner 
boundary.  Important  characteristics  in  grid  generation  are 
the  ability  to  specify  the  spacing  between  mesh  points  at  the 
boundary,  in  the  direction  normal  to  the  boundary,  and  the 
control  of  the  angles  with  which  mesh  lines  intersect  the 
boundaries  which  is  known  as  orthogonality. 
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The  grid  size  for  the  sonic  surface  (M^l.O)  is  107  x  60  and 
for  the  supersonic  surfaces  (M^l.l,  1.2)  is  115  x  60  as 
shown  in  Fig. (18) .  Since  the  surfaces  are  symmetrical,  half  of 
the  grid  or  eventually  the  lower  surface  was  rotated  for  11 
planes  plus  two  more  for  symmetry  to  generate  an  axi- 
symmetric  after  body  surface. 

The  OVERFLOW  program  [Ref  .10]  developed  by  Ames  Research 
Center  which  uses  3-D  Navier- Stokes  and  Euler  solver  for 
viscous/  inviscid  flow  was  applied.  The  results  are  shown  in 
Figs.  19-21  for  three  different  boundary  surfaces  at  M^  =  1.0, 
1.1,  and  1.2  with  1500  iterations.  It  was  found  that  the 
density  residuals  decreased  by  more  than  one  order  of 
magnitude  over  1500  iterations.  After  3000  iterations,  the 
solution  converged  by  two  orders  of  magnitude. 

In  case  of  the  sonic  surface  Fig.  (19)  with  M^l.0,  the 
Mach  lines  contour  shows  the  shock  is  forming  downstream  at  a 
local  Mach  M=1.65  at  which  the  flow  becomes  subsonic 
downstream  of  the  shock.  This  result  complies  with  the  small 
perturbation  transonic  solution  obtained  earlier  in  Fig.  (15)  . 
In  other  words,  the  flow  is  shock  free  over  the  sonic  surface 
in  for  M^l.0. 

For  the  supersonic  boundary  surfaces  NL^l.10,  and  1.20,  we 
can  see  that  the  flow  starts  at  M=1.10  and  M=1.20  respectively 
and  follows  the  afterbody  surface  until  it  shocks.  The  shock 
is  formed  at  local  Mach  M=1.75  and  1.85  respectively  as  shown 
in  Figs.  20  and  21. 
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On  the  other  hand,  these  results  agree  with  the  small 
perturbation  solution  over  the  transonic  range  as  plotted  in 
Figs.  11  and  13  where  the  local  Mach  is  represented  by  the 
dash  line  and  it  is  increasing  towards  the  afterbody  surface 
until  it  shocks.  This  indicates  that  the  flow  is  shockless 
over  these  supersonic  boundary  surfaces  in  the  transonic 
regime . 

The  velocity  vectors  are  plotted  in  Fig.  (22)  for  M^l.O 
which  shows  that  the  flow  is  following  the  afterbody  surface 
until  the  shock,  and  then  flow  separation  takes  place  due  to 
the  steep  pressure  rise  and  the  introduction  of  numerical 
viscosity  into  the  flowfield  at  this  location. 

These  boundary  surfaces  obtained  are  of  interest  for 
design  of  body  of  revolution  such  as  missiles  afterbody  or 
aircraft  engines,  in  which  a  patching  technique  may  be  applied 
to  these  surfaces  to  get  the  best  shockless  surface  in  the 
transonic  range. 
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VII.    CONCLUSION  AND  RECOMMENDATION 

We  have  shown  in  this  research  a  solution  to  the  non- 
linear transonic  small  perturbation  equation  by  implementing 
the  separation  of  variables  technique.  The  x- dependence  was 
integrated  numerically  and  the  r- dependence  was  solved  by  an 
outer  expansion  method.  It  was  found  that  the  parameter  (1- 
mJ2)  has  strong  effect  in  specifying  the  shape  of  the  boundary 
surfaces  of  interest  for  design.  A  3-D  viscous/  inviscid 
Navier- Stokes/  Euler  solver  in  Computational  Fluid  Dynamics 
(CFD)  confirmed  our  results  obtained  from  the  small 
perturbation  technique  for  the  boundary  surfaces  at  M^l.O, 
1.10,  and  1.2  0  .  In  other  words,  we  achieved  our  goal  of 
designing  shockless  surface  in  the  transonic  range. 

This  research  may  be  extended  using  a  patching  technique 
to  obtain  the  best  shockless  surface  in  the  transonic  flows. 
Patching  criteria  will  have  to  be  developed  based  on  the 
mathematics  and  the  flow  constraints.  Also,  connection  between 
small  perturbation  solution  and  CFD  might  eventually  be  done 
internally  in  the  computer.  The  small  perturbation  is  the  most 
efficient  method  to  define  or  to  design  axi-symmetric 
afterbody  transonic  surfaces,  however  CFD  can  be  used  to 
predict  the  performance  of  these  aerodynamic  surfaces. 
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APPENDIX  B 
COMPUTER   PROGRAMS 


**************************************************************** 

*  * 

*  This  program  is  designed  to  calculate   KSI(X),  * 

*  K(X),  and  X  using  numerical  integration  based  * 

*  on  trapezoidal  rule  to  solve  Eqn. (24)  and  * 

*  plotting  the  output  as   shown  in  Figs . 1  &  3  .  * 

*  * 
***************************************************** 

PROGRAM  KSI 

REAL  M (3) ,X (0:401, 3) , A, Al , B (401 , 3 ) 
+      ,KSI,P,FUNC,H,DD,Q,L,K(401,3) ,N,N1 
INTEGER  YY 
OPEN (UNIT=9,FILE=' KSI' , STATUS =' UNKNOWN' ) 

PRINT  *, 'ENTER  LOWER  &  UPPER  BOUND,  #  OF  INTERVALS, 
+#OF  DATA  SET' 
READ  *,  A,  Al,  N,  Nl 

100    DO  10  1=1,3 

PRINT  * , ' ENTER  MACH  NO .  ' 

READ  *,  M(I) 
C       M(I)  =  0.9+1*0.1 

IF  (M(I) .LT.1.0)  THEN 

P=-l. 

ELSE 

P=l. 

ENDIF 

DO  2  0  J=1,N1 

B(J,I)  =  J*A1/N1 

H  =  (B(J,I) -A)/N 

AREA  =  0. 

K(J,I)=(-0.0102/(M(I)**2) )* ( (ABS( (B(J,I)**2) - 
+  (  (ABS(1-M(I) **2)  ) ** 1.7574)  )  ) **  (2./3.)  + 
+ (ABS(1-M(I) **2))**1.1716) 

DO  30  L=l  ,  N 

KSI  =  A  +  (L-0.5) *H 

DD  =  (KSI**2+P* (ABS(1-M(I) **2) )**1.7574) 
IF  (DD.GT.0.)  THEN 
FUNC  =1./ (DD** (1./3. ) ) 

Q  =  1. 
ELSE 
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FUNC  =  l./( (-l*DD)**(l./3.) 

Q  =  -1. 

ENDIF 

AREA  =  AREA  +  H*FUNC*Q 

X(J,I)  =  AREA 

30 

CONTINUE 

20 

CONTINUE 

10 

CONTINUE 

DO  40  J=1,N1 

PRINT  50,  (B(J,I) ,X(J,I) ,K(J,I) ,1=1,3) 
50     F0RMAT(1X,3 (1F6 . 3 , 2F10 . 5 , IX) ) 
40     CONTINUE 

200    PRINT*, 'TRY  AGAIN  ?  ENTER  1  FOR  YES , 2  FOR  DATA 
+FILE, OTHERS  FOR  NO' 
READ  *,  YY 
IF  (YY.EQ.l)  THEN 

GO  TO  100 
ELSE 

IF  (YY.EQ.2)  THEN 
DO  110  J=1,N1 

WRITE  (9,50)  (B(J,I) ,X(J,I) ,K(J,I) ,1=1,3) 
110  CONTINUE 

GO  TO  200 
ENDIF 
ENDIF 
END 
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**************************************************************** 

*  * 

*  This  program  is  written  to  calculate   ZETA(r)  and       * 

*  r  using  Eqn. (25)  and  plotting  the  output  in  Fig. 3  .      * 

*  * 
**************************************************************** 


program  zeta 

real  M(3) , zeta (0 : 14 , 3 ) , r (14) 

open (unit=9 , f ile=' zz' ; status=' unknown' ) 

do  10  1=1,3 
c       M(I)=0. 7+1*0.1 

print  *,  'enter  mach  no.' 
read  *,  M(I) 

do  20  J=l,14 

zeta(0, I) =1000 

print  *, 'enter  r  values' 

read  *,  r(J) 
c       r(J)=0.l*J 

zeta (J, I)  =  (4/ (r(J) **2)  )+  (abs(l- (M(I) **2)  )  * 
+ (r(J) **2.8284) ) + ( (1- (M(I) **2) ) **2) * (r(J) **7.657) /50.63 

if  (zeta(J, I) .GE.zeta(J-l, I) )  zeta ( J, I) =zeta ( J-l, I) 
20      continue 
10     continue 

print  30,  M 
30      format (llx,5F10. 4) 

do  40  J=l, 14 
c       print  50,  r ( J) , (zeta ( J, I) , 1=1 , 4) 

write  (9,50)   (r ( J) , zeta ( J, I) , 1=1, 3 ) 
50      format  (4x, 3 (2F10 . 4 , 2x) ) 
40      continue 

end 
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**************************************************************** 

*  * 

*  This  program  is  written  to  calculate   Z(r)  and   r        * 

*  using  Eqn. (59)  and  showing  the  output  in  Fig. 4  .         * 

*  * 
**************************************************** 


PROGRAM   ZR 

REAL  M(4) ,R(1000) ,Z(0:1000,4) 

0PEN(UNIT=99,FILE='ZR1' , STATUS =' UNKNOWN' ) 

DO  10  1=1,4 
C       M(I)  =  0.7  +  1*0.1 

PRINT  *,  'ENTER  MACH  NO.  ' 
READ  *,  M(I) 

DO  20  J=l,1000 
Z(0,I)  =  0.0 
R(J)  =  0.001*J 

Z  (  J  ,  I  )     =     (R(J)**3)/(8 

(2.8284* (R(J) **4.82  84) *ABS ( 1-  (M ( I) **2  )  )  )  - 

+(0.1512* (R(J) **9.657) * ( (1-M (I) **2 ) **2) ) ) 
IF  (Z(J,I) .LE.0.0)   Z(J,I)=Z(J-1,I) 
20     CONTINUE 
10     CONTINUE 


DO  40  J=l,1000 
C       PRINT  50,  R(J)  ,  (Z(J,I) ,1=1,4) 

WRITE  (99,50)  R(J),  (Z(J,I),  1=1,4) 
50     FORMAT(4X,5F10.6) 
40     CONTINUE 

END 
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**************************************************************** 

*  * 

*  This  program  is  written  to  determine  Rmin  for  * 

*  transonic  Mach  numbers  using  Eqn. (60)  and  the  * 

*  output  is  shown  in  Fig. 5  .  * 

*  * 
**************************************************************** 


PROGRAM  RMIN 

REAL  M,R 

OPEN (UNIT=88,FILE=' RMIN' , STATUS =' unknown' ) 

DO  100  M  =  0.8,1.2,0.0001 
IF(M.LT.l.O)   THEN 

R=1.2074/ ((1-M**2)**0.2071) 
ELSE 

IF(M.GT.l.O)   THEN 

R=1.2074/  (  (M**2-l) ** 0.2 071) 
ELSE 

R  =  7.0 
ENDIF 
ENDIF 
C       PRINT  80,  M,R 

WRITE  (88,80)  M,R 
80     FORMAT(1X,F7.5,7X,1F10.6) 
100    CONTINUE 
END 
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**************************************************************** 

*  * 

*  This  program  is  designed  to  numerically  integrate        * 

*  Z(r)  ,Eqn.(64),  based  on  trapezoidal  rule  .  * 

*  * 
**************************************************************** 


PROGRAM   IZR 

REAL  M(3) ,R(401,3) , A, Al , IZ (0 : 401 , 3 ) , H, FUNC, N, Nl , Ro 

INTEGER  YY 

OPEN(UNIT=55,FILE='IZR' , STATUS =' UNKNOWN' ) 

100    DO  10  1-1,3 

C       M(I)  =  0.7  +  1*0.1 

PRINT  *,  'ENTER   MACH   NO.  ' 

READ  *,  M(I) 

PRINT  *,  'DETERMINE  THE  VALUES  OF  MACH  #',M(I) 


PRINT  * 
READ  *, 
PRINT  * 
READ  *, 
PRINT  * 
READ  *,N 
PRINT  * 
READ  *, 
PRINT  * 


'ENTER  Rmin   ,A' 
A 

'ENTER  THE  UPPER  LIMIT  OF  R  ,A1' 
Al 

'ENTER  #  OF  INTERVALS  ,  N' 

'HOW  MANY  DATA  SET  DERIVED  ?  Nl ' 
Nl 
'INPUT  Ro  FOR  THIS  MACH  #  ,  Ro' 


READ  *,  Ro 

DO  20  J=1,N1 
R(J,I)  =  J*A1/N1 
H  =  (R(J,I) -A)/N 
AREA  =  0. 

DO  30  L=1,N 
RR  =  A+ (L-0.5) *H 

FUNC= (RR**3) / (8- (2.82  84* (RR**4.82  84) *ABS (1- (M (I) **2) ) ) 
+  - (0.1512* (RR**9.657) * ( (1-M (I) **2 . ) **2 . ) ) ) 
AREA  =  AREA  +  H*FUNC 
IZ(J,I)  =  Ro  -  AREA 


30 
20 
10 


CONTINUE 
CONTINUE 
CONTINUE 


DO  40  J=1,N1 
PRINT  50,  (R(J,I) ,IZ(J,I) ,1-1,3) 
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WRITE (55, 50)   (R ( J, I) , IZ ( J, I) , 1=1 , 3 ) 
50     FORMAT  ( IX, 3 (F6 . 3 , F9 . 6 , 3X) ) 
4  0     CONTINUE 

PRINT  *,'TRY  AGAIN  ?  ENTER  1  FOR  #YES#,  OTHERS  FOR  #NO# ' 

READ  *,  YY 

IF  (YY.EQ.l)  GO  TO  100 

END 
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**************************************************************** 


This  program  is  using  numerical  integration  to 
determine  the  supersonic  boundary  surfaces  , 
Eqn.  (64) ,  and  calculating  pressure  coefficient 
and  local  Mach  number  profiles,  using  Eqns.(40 
&  49) .  The  output  is  shown  in  Figs.  (6-13)  . 


**************************************************************** 


PROGRAM  CPP 

REAL  X(401, 9) ,Y(401,6) 

INTEGER  I#JfK,L, IWRITE 


XR(401) ,YR(401) ,M,MX,XYMAT 


OPEN (UNIT=52 , FILE= ' cpl5 . dat ' , STATUS= 
OPEN (UNIT=53 , FILE= ' cpll . dat ' , STATUS= 
OPEN (UNIT=54 , FILE= ' cpl2 . dat ' , STATUS= 
OPEN (UNIT=77 , FILE= ' cpp . dat ' , STATUS= ' 
OPEN (UNIT=78 , FILE= ' ml05 . dat ' , STATUS= 
OPEN (UNIT=79 , FILE= ' tirel . dat ' , STATUS 
OPEN (UNIT=80 , FILE= ' tire2 . dat ' , STATUS 


'UNKNOWN' ) 
'UNKNOWN' ) 
'UNKNOWN' ) 
UNKNOWN' ) 
'UNKNOWN' ) 
=' UNKNOWN' ) 
=' UNKNOWN' ) 


100    OPEN (UNIT=1,FILE='KSI. DAT' , STATUS= ' OLD' ) 
OPEN (UNIT=2,FILE='IZR. DAT' , STATUS= ' OLD' ) 
REWIND (UNIT=1) 
REWIND (UNIT=2) 
IWRITE  =78 
MACP=52 


READ  (1,*,END=50) 
50  READ  (2,*,END=60) 
60     DO  10  J=2,6,2 


( (X(I,J) ,J=1,9) ,1=1,401) 
( (Y(K,L) ,L=1,6) ,K=1,401) 


PRINT  *,'  Ro  is  rmin  ;  EPSILON   is  the  accuracy 

PRINT  *, 'ENTER   MACH   NO.,  Ro  ,  EPSILON' 

READ   *,  M,RO,EPS 

WRITE  (77,*)  M,RO,EPS 

XM  =  1-M**2 

L  =  J/2*3 

DO  30  I  =  1,401 
DO  20  K  =  1,401 

IF(Y(I, J) .LT.0.001)  THEN 

GO  TO  10 
ELSE 

XYMAT  =  ABS(Y(I, J) -ABS(X(K,L) ) ) 
IF (XYMAT.GT.EPS)  THEN 
GO  TO  20 
ELSE 
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XR(K)  =  X(K,L-1) /RO 
YR(I)  =  Y(I, J-l) /RO 

ZETA  =  (4./ (Y(I, J-l)**2) )  +  ABS(XM) *Y(I, J-l) **2. 8284 
+  +  XM**2* (Y(I, J-l) **7.657)/50.63 

Cp  =  (-2/ (M**2*2.4) ) * (0.2208*ZETA* (X(K,L-2) **2  + 
+  (ABS(XM) )**1.7574)**(l./3.)+XM) 

MX  =  SQRT(ABS( (M**2* (1-Cp)/ (1+0 . 2*M**2*Cp) ) ) ) 
C  PRINT  90,  X(K,L-1) ,Y(I, J-l) ,Cp,MX 

WRITE (77, 40)  X(K,L-1) ,Y(I,J-l),Cp,MX 
WRITE (MACP, 40)  XR(K)  ,YR(I)  , Cp  ,MX 
WRITE (IWRITE, 70)  XR(K)  ,  YR(I) 
GO  TO  3  0 
ENDIF 
ENDIF 

20     CONTINUE 
3  0     CONTINUE 

IWRITE=IWRITE+1 

MACP  =MACP+1 
10     CONTINUE 

40     FORMAT  ( IX, 4 (F15 . 7 , 2X) ) 

70     FORMAT  (1X,2F17.9) 

90     FORMAT  (IX, 4 (F15 . 7 , 2X) ) 

PRINT  * , ' TRY  AGAIN  ?  1  FOR  YES ! , 2  FOR  NO ! ! ' 
READ  *,  IY 

IF(IY.EQ.l)  GO  TO  100 
END 
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**************************************************************** 

*  * 

*  This  program  is  written  to  determine  the  sonic  * 

*  boundary   surface  M  =  1.0,  ro  =  4.0  using  Eqn. (78) ,       * 

*  and  calculating  pressure  coefficient  and  local  * 

*  Mach  profile  using  Eqns . (79  &  80) .  The  output  is         * 

*  shown  in  Figs.  16  &  17  .  Some  output  files  are  used      * 

*  in  CFD.  * 

*  * 
**************************************************************** 


program  sonic 

real  M(0:160) ,r (0:160) ,x (0:160) ,  zeta (0:160) ,ksi (0:160) , 
+      z(0:160)  ,k(0:160)  ,Cp(0:160)  ,    RR (0 : 160) , YY (0 : 160) , 
+      xx(0:160) ,xxx(0:160) 

open (unit=ll, f ile=' sonicl .dat ' , status=' unknown' ) 
open (unit=22 , f ile=' sonic2 .dat ' , status=' unknown' ) 
open (unit=33 , f ile=' sonic3 .dat' , status=' unknown' ) 
open (unit=66 , f ile= ' sonic4 . dat ' , status= ' unknown' ) 

c  sonicl.dat  =>  x  ,  r  ,  Cp  ,  M 

c  sonic2.dat  =>  x  ,  r   (horintal) 

c  sonic3.dat  =>  x  ,  r   (vertical) 

c  sonic4.dat  =>  x/ro  ,  r/ro  ,  Cp  ,  M 

do  10  I  =  0,  159 

x(I)  =  1*0.1 

ro  =  4.0 

r(I)  =  (abs(ro**4.  -  0 . 004* (x(I) **4 . ) ) ) **0  .25 

RR(I)  =  r(I)/ro 

YY(I)  =  -  1.0*  RR(I) 

xx(I)  =  x(I)/ro 

xxx(I)  =  3.975  -  xx(I) 

zeta(I)  =  4./ (r(I) **2. ) 

ksi(I)  =  (x(I)**3.)/27.0 

z(I)  =  r(I) **3. 

k(I)  =  -  0.00101* (x(I) **4.) 

Cp(I)  =  -  0.1840  *  zeta(I) * (ksi(I) **0.6667) 

M(I)  =  (abs( (l-Cp(I) )/(l  +  0.2*Cp(I) ) ) )**0.50 
c      print  *,  xx(I)  ;RR(I)  ,Cp(I)  ,M(I) 
10    continue 

do  20  I  =0,  159 

write(ll,50)  xx(I)  ,  rr(I)  ,  Cp(I)  ,  M(I) 
write(22,60)  xxx(I)  ,  RR(I)  ,  YY(I) 
write(33,70)  xxx(I),YY(I) 
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write(66,5C)  x(I)  ,r(I)  ,Cp(I)  ,M(I 
50     format (2x, If 8.4, 2x,3f 14.7, 2x) 
60     format (2x, If 8.3, 2x,2f 12.6, 2x) 
70     format  (lOx,  If  12  .  6  ,  5x,  f  12  .  6) 
20     continue 

write (34, 80) (xxx(I) ,i=0,159,3) 
write (34, 80) (xxx(i) ,1-156,0,-3) 
write(34,80) (yy(i) ,i=0,159,3) 
write (34, 80) (rr (i) , i=156,0, -3) 
80     format (5 (lx,f 12.6, ',')  ) 

do  30  I  =158,  0  , -1 
write(33,70)  xxx(I),RR(I) 
30     continue 
end 
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